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THEORET I CAL.. "WORK 

’ ** fj 

Personnel 

R. W. Thresher has completed his work on the contract and has received 
his Ph.D. Wo student is presently employed on the theoretical portion of 
the work. 

Progress 

The computer program has been run several times for the case of uni- 
form loading on the crack shown in Figure 1, Stress intensity factors have 
been computed for values of D ranging from D~0.0toD=0.7 and for values 
of A/1 ranging from A/1 = 0.0 to A/1 = 0.85. 

The results of these computations are presented in Dr. Thresher's 
thesis, a copy of which is attached. Figures 6, 12, 13 and 14 present 
stress intensity as a function of position along the crack border. Figure 15 
presents the back surface magnification ' factor as a function of A/1. Com- 
parisons of theoretical and experimental work are shown in Figures 16 and 


EXPERIMENTAL WORK 


Personnel 

A master’s degree candidate has been hired 1/2 time to work on this 
portion of the research. 

Progress 

Testing has continued on the width effect. In excess of 60 data points 
have been added to the data presented in the last report. Figures 2 through 



11 show all the width effect data taken to date. The hollow points show 
the most recent data while the solid points show the data presented in 
the last report. 

The apparent K.^ (AK^) in Figures 2 through 11 is still defined as 



Now that enough data are avilable, a new apparent K IC will be defined, 
namely 

. - wfif 

as was discussed in the last report. This new apparent K^ c will he used 
to compute width effect correction factors. 

In the last report- it was mentioned that the width effect may have 
influenced some of the depth effect data. It appears now that this is 
not the case since the depth effect specimens were purposely made with 
large values of W/2C. The following criterion were applied to the depth 
effect data: 

invalid if 

W/2C <4.0 and a/2C > .30 
W/2C <2.0 and a/2C < .30 

This resulted in eliminating only four of the total depth effect data 


points. 
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ABSTRACT 


A SURFACE CRACK IN A FINITE SOLID 

A solution to the problem of a circular crack partially 
embedded in a solid of finite thickness is presented, A super- 
position and iteration technique is used to deter mine the stress 
intensity factor numerically. The stress intensity factor is 
determined as a function of position around the crack front for 
a variety of crack depths and thicknesses. The results of this 
study are compared to experimental data for a semi-elliptical 
surface flaw in a brittle material. 

In addition, a solution for the partially closed Griffith 
crack is presented. This formulation leads to a closed form 
solution for the stress and displacement on the plane of the crack. 

A method for finding the open crack length is presented. The 
problem of pure bending which closes one end of the crack is solved 
and the stress intensity factor is calculated. 

Robert Wallace Thresher 
Department of Mechanical Engineering 
Colorado State University 
Fort Collins, Colorado 80521 
August, 1970 
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INTRODUCTION TO THE THESIS 


The present form of fracture mechanics started with the work 
of Griffith in 1920 (1). The essence of Griffith’s work stated that 
a crack will propagate if by doing so it can lower the total energy 
of the system. The stress analysis used by Griffith to caleulate the 
stored elastic energy was taken from the work of Inglis (2) for an 
elliptic hole in the center of an infinite plate, where the minor 
axis of the ellipse was allowed to approach zero. 

Another step toward the present theory of fracture mechanics 
was made when Irwin (3) pointed out that for metals the Griffith-type 
energy balance must be between the stored strain energy and the surface 
energy plus the work done during plastic deformation. Irwin also 
noted that for most engineering materials the work done against sur- 
face tension is small compared to the plastic work. 

The Irwin-Grif fith energy balance for a crack would appear as 
shown in Figure 1. As the crack size increases the stored elastic 
energy decreases and the amount of plastic work necessary to propagate 
the crack increases. However, the stored strain energy decreases at 
a faster rate and therefore creates an unstable situation when the 
crack is larger than a certain critical size. 

In 1957 Irwin (4) showed that the strain energy release rate 
is related to the stress intensity factor K^. where the stress inten- 
sity factor can be defined as 


1 
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K t = Limit /2 tt6 a 
-*• * . n XX „ 

o -*• 0 x=0 

for the situation illustrated in the bottom half of Figure 1. Since 
the strain energy release rate and the stress intensity factor are • 
related, it is permissible to talk about a critical stress intensity 
factor. The critical stress intensity factor is the value of the 
stress intensity factor at which unstable crack propagation occurs. 

The critical stress intensity factor, which is called the fracture 
toughness, has been shown by experiment (5) to be a material property. 
The fracture mechanics failure criterion states that a crack will 
propagate if K-j- > Kjq , where K is the fracture toughness. 

To determine whether a cracked component will fail in a given 
loading situation, one must have the fracture toughness for the 
material in question, and an elastic stress analysis from which to 
calculate . For example, consider a two-dimensional crack of 
length 2a centered in a large plate loaded with a uniform tension 
a . The plate is a brittle material which has a fracture toughness of 
550 psi /Ln . A stress analysis is performed to obtain a ^ on x = 0, 
from which the stress intensity factor is calculated. This computation 
gives 

Kj = a/rra 

The fracture mechanics failure criterion states that this plate will 
fail when 

a/rra > 550 

This equation gives the maximum allowable load O for a given crack 
length 2a . For most engineering materials the fracture toughness is 
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well documented. The difficult part of the problem Is the stress 
analysis. Solutions to some of the more common problems can be found 
in the literature (5). However, many important problems are still 
unsolved. For example, surface cracks are very often seen in structural 
components, but only limited analysis has been performed on this type 
of crack. 

The purpose of this thesis is to solve two problems in linear 
elastic fracture mechanics to obtain the stress intensity factor. 

Part I presents the solution for a surface flaw in a finite thickness 
plate. Part II presents the solution for a two dimensional through 
crack in which the applied load closes one end of the crack. 



INTRODUCTION TO PART I 


The three dimensional analysis of cracks started in 1945 with 
the work of Sneddon (6) , who solved the problem of the circular crack 
in an infinite solid, where the crack was opened by a constant pres- 
sure. The corresponding problem for the elliptical crack in an 
infinite solid was done by Green and Sneddon in 1960 (7). More recently, 
Kassir and Sih (8, 9) solved the embedded elliptical crack for a pre- 
scribed shear, and later for a linearly varying normal pressure. The 
problems of an embedded hyperbola and an embedded parabola were solved 
by Shah (10) . In addition, Shah solved the embedded elliptical crack 
for a normal pressure loading specified by a restricted polynomial (11) . 

Only limited analytical work has been done, on surface crack', 
problems. In 1962 Irwin (12) using the solution of Green and Sneddon 
estimated the stress intensity factor for a semi-elliptical crack in 
the surface of a half space. In 1965 Smith (13) solved the semi- 
circular crack in the surface of a half space, where the crack surface 
was loaded with an arbitrary normal pressure. In a continuation of 
Smith’s work. Smith and Alavi solved the problem of a circular crack 
embedded in a half space (14) and the problem of the circular crack 
only partially embedded (15, 16). 

The work contained herein is an extension and a refinement of 
the work of Smith and Alavi. It is the purpose of this portion of the 
thesis to present the solution to the part-circular crack in a finite 
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thickness solid. The work centers around determining stress intensity 
factors for circular flaws embedded an arbitrary depth into a plate 
of finite thickness. 



CHAPTER I 


THE SURFACE FLAW PROBLEM 


1.1 The Problem Statement 

Consider a circular crack embedded only part way into a plate 
as illustrated in Figure 2. The problem is to determine the stress 
intensity factor as a function of 0 for a variety of depths D and 
A/T ratios under uniform tension loading. 

Again referring to Figure 2, the boundary conditions may be 
stated as follows: 

1. All normal and shearing stresses must vanish on both the 
front and back surfaces. 


2. On the plane z=0: 

(a) a » - P(r,0) 

ZZ 2=0 

■'(b) u =0 

2=0 


(c) 




0 


inside the crack, 
outside the crack, 
all 2=0. 


The solution to this problem will be superimposed with a solution 
having a normal stress + P(r,0) inside the crack to arrive at the 
boundary condition a =0 inside the crack. 

* 77 

ZZ 2=0 


6 
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1.2 The Solution Method 

The solution method which will be used is the Schwarz alter- 
nating technique (17). In this particular case the technique consists 
of using two elastic solutions coupled with an iteration procedure 
programmed for the digital computer to solve numerically for the stress 
intensity factor as a function of position along the crack border. 

The solutions used in this procedure are: 

1. The circular crack in an infinite solid where the crack 
surface pressure is prescribed in the form of a Fourier 
series; 

2. Stress in a half space due to pressure and shear on a 
rectangular portion of the half space boundary. 

The iteration procedure works as follows: 

1. The crack pressure, prescribed in the form of a Fourier 
•series, is applied to the crack surface. For the case of 

simple tension, a constant pressure is applied to the 
crack surface. Solution 1 is used to calculate the result- 
ing stress in the solid at the desired location for the 
front surface. 

2. Stresses at the location of the front surface are removed 
using solution 2. This is accomplished by dividing the 
front surface into many small rectangles and applying 
"freeing" stresses to each rectangle which are equal and 
opposite to the stresses at the center of the rectangle as 
computed by solution 1. This freeing of the front surface 
produces a residual stress on the crack surface, which is 
computed using solution 2. 
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3. The negative of the stress residual on the crack surface 
is approximated with a Fourier series and this loading is 
then applied to free the crack surface of stress. Solution 
1 is again used to compute the stress on the front surface. 

4. Steps 2 and 3 are repeated until the residuals become 
negligible. 

As the iteration continues, the stresses due to each step are 
also computed for the back surface and collected in a running sum. 

After convergence with the front surface has been obtained the iter- 
ation is done between the crack and the back surface, while a running 
total of residuals is kept on surface 1. This whole procedure is then 
repeated until, residuals are negligible. 

The elastic solution for the circular crack in an infinite 
solid will be developed in Chapter II and the solution for a rectangle 
on the surface of a half space will be discussed in Chapter III. 



CHAPTER II 


THE CIRCULAR CRACK IN AN INFINITE SOLID 
SOLUTION 1 


2.1 Boundary Conditions 

Consider an elastic solid with a circular crack of radius "a" 
located on the plane z=0 which is opened by normal pressure P(r,0), 
The crack and coordinate system are shown in Figure 3. The boundary 
conditions for this problem may be written as follows: 


(a) 

cr 

zz 

= - P(r,8) 

1 3=0 

for 0 < r < a 

(2.1.1a) 

(b) 

U 

= 0 


for r > a 

(2.1.1b) 


z 

2=0 




(c) 

T 

rz 

z=0 = T0Z 

= 0 
z=0 


(2.1.1c) 

(d) 

a. ■ 
ij 

-»• 0 as x,y -*■ “ 


(2.1. Id) 


In addition, only normal crack pressures symmetric with the plane 
z=0 will be considered. 

2.2 The Potential Formulation 

The three dimensional form of the Navier equations can be 

written 

G V 2 U^* (A + G) U ± ± = 0 (2.2.1) 


9 
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For the special class of problems where the shearing stresses 
are zero on the plane z=0 the following potential formulation satisfies 
the Navier equations: 


where 


+ 
9r ' 

(f5 2 - 

l)z 

3 2 (J) 

9r9z 


l u 

r 90 

+ (8 2 

- 1) 

z 

r 

9 2 <f> 

909z 

- e 2 

M + 

9z + 

(8 2 - 

l)z 

9 2 <j> 

9z 2 

2(1 - 

- v) 




(1 - 

2v) 





( 2 . 2 . 2 ) 


provided that V 2 4> = 0 . This formation is presented by Green 

and Zerna (18) . 

Using Hooke's Law the stress components in terms of the 
potential function are easily found as 


a 


rr 


cr 


09 


a 

zz 


- 2 (f3 2 -2) + 2 ^ + 2(8 2 -1)z r-^- 

9z 2 9r 2 9r 2 9z 


- 2 (8 z -2) ^L + 2 |£ + -L 


9z 2 r 3r r 2 90 2 


+ 2(8 2 -l) 2 (8 2 -l) 

r 2 90 2 9 z r 3r3z 


- 2(8 2 -l) 2(6 2 -1 )z 


9z^ 


9z c 
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re 


2 3 2 4 > 
r 3r30 


+ 2(3 2 -l) - 


3 3 4> 

303z3r 


_ 2 _ 3 £ 

2 90 


- 2(0 2 -l) 


3 2 cJ) 

303 2 


T 

rz 


2(3 Z -l)z 

3r3z 2 


T = 2 (3 2 “1) f (2.2.3) 

02 r 303z 2 

The displacements in (2.2.2) have been nondimensionalized through 
division by "a", the crack radius. The stress has been nondimension- 
alized with respect to the shear modulus. 

From the form of the potential formulation it can be seen that 
the boundary condition (2.1.1c) is automatically satisfied while 
boundary conditions (2.1.1a) and (2.1.1b) on the plane z=0 reduce to 


3 z <fr = P(r,9) 

3z 2 2(3 2 - 1)' 


3i 

3z 


0 


for 0 < r < a 


for r > 0 


(2.2.4a) 


(2.2.4b) 


2.3 The Potential Function 

To satisfy the mixed boundary conditions of (2.2.4a) and 
(2.2.4b) the potential function, <J> , will be assumed in the form 


N 

= l cos 
n=0 



J n (£r) e“ z? d£ 




(2.3.1) 
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This potential function clearly satisfies the boundary condition (2.1. Id), 

_ z r 

since both J (£r) and e tend toward zero for large r and z . 
n 

The remaining boundary conditions to be satisfied are (2.2.4a) and 
(2.2.4b), which after substl t-ni-irm for $ may be rewritten as 


P(r,6) 

2 ($ 2 -l) 


0 


N 

£ cos (n0) 
n-0 


00 

I 5f(5) j 

J n n 


(Cr) d£ (2.2.3a) 

for r < 1 


00 

N r 

I COS (n8) / f n (C) J n (Cr) d£ 
n=0 J * 


(2.3.2b) 


for r > 1 


To find the unknown function f n (£), the funccxon P(r,0)/2(g 2 -l) 
will be expanded in a Fourier cosine series as follows: 


P(r.,8) 

2(g z -l) 


N 

\ B (r) cos (n8) 
n=0 n 


(2.3.3) 


where the Fourier coefficients B Q (r) are given by 


TT 

B 0 (r) = ~ f P(r,0) d0 

2tt(8 2 -1) j q 

and- ^ (2.3.4) 

B (r) = f P(r,0) cos (n0) d0 

7r(3 2 -l) J q 


Now the boundary conditions (2.3.2a) and (2.3.2b) may be further 
simplified to give the following set of dual integral equations: 
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CO 

5 n (r) = f ?f n (?) J n (?r) d? 


(2.3.5a) 


for 0 < r < 1 


OO 

o = j f n (D J n (?r) d£ 


(2.3.5b) 


for r > 1 


The solution to these dual Integral equations for the unknown function 
f n (£) , due to Busbridge, is given by Sneddon in (19) as 


1 l 

E »® - \[W f W 5 ' 1,dri / B n 

» *n n 


n+1 

(tip) dp 


(2.3.6) 


2.4 The Function f (£) 

n ^ 


In order to reduce the complexity of the integrations required 
to determine f^(£) , the function B^(r) will be expanded in the 
power series 


B (r) = 7 C p 

n -n n 

P 0 


(2.4.1) 


where C = constant. Substitution of this representation for B^(r) 
into (2.3.6) and performing the integration on p gives 




l 

" ^ jo ? " P H " P / r,P+% w e ' )4n - <2 - 4 - 2) 


where 



1 fp +n+ 3 ~ 

i .2 . 


(2.4.3) 
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To this point the solution development has followed the work of Smith 
and Alavi (15) , who expanded the Bessel function in a power 

series which was then integrated term by term and substituted into 
the potential function $ to obtain a series solution for the stress 
components. Because the series solution did not converge for 
r + z < 1 , they were forced to developoa closed form solution for 
small r and z . The main drawback with the closed form solution 
was that it required a considerable amount of hand computation to 
complete the required integrations. Even then it took a significant 
amount of computer time to calculate numerical values, for the stresses. 
As a result of these difficulties, Alavi' s solution is limited to 
three terms in both the Fourier series (2.3.3) and the power series 
(2.4.1). In order to extend the solution to include more terms in the 
Fourier series and power series, a numerical integration scheme was 
devised. 

Watson (20) gives the integral representation for J I1+ ^(£ T l) as 

J . ,(5n) = -&&>!!!? f dt (2 . 4 . 4) 

* r(n+i)/fT 

Substitution of this integral representation into equation (2.4.2), and 
then substitution of equation (2.4.2) into equation (2.3.1) gives the 
result 

= ? (-8) y c P H P . X 

n£o TT2 n+1 r(n+l) p-0 11 n 

11 00 

1 1 riP+a+2<1 " t2)n<itdr| f S n J n (?r)e _ ^ (z “ iT1t:) d£ 
*0-1. 0 (2.4.5) 
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The last integral in (2.4.5) has the value 


OO 




r(2n+l)r 


n 


n , , (2n±l)/2 


2 f (n+l) [ (z-int) +r ] 


(2.4.6) 


Substitution of equation (2.4.6) into equation (2.4.5) gives 


N 


l ±1 ^ J c'.'r'K.,,) 


n=0 /F 2 2n+1 |T(n+l)] 2 p=0 n • n 


(2.4.7) 


where 


I(n,p) 


1 1 

If 

0-1 


r^p+n+2 a ^n 
t(z-iht) 2 +r 2 ] (2n+1)/2 


dtdn (2.4.8) 


2.5 The Stress Components 

To determine the stress components, the equation (2.4.8) for 
I(n,p) must be substituted into the equations (2.2.3), which gives 
the following integrals to be evaluated: 


FI(1) 


1 1 


■ff^d 


(2n+l)/2 


dtdn 


0 -1 


1 1 


FI (2) = // (2n+3)/2 dtdf| 


Oi' -1 


1 1 


FI(3) = // (2n+5)/2 dtdn 


0 — 1 


(2.5.1) 


( 2 . 5 . 2 ) 


(2.5.3) 
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where 


and 


fi( 4> - // dtan 

0 -1 4 

FI(5) - // TtfsSfk dtdn 
0 -1 4 

FI(6) ’ // jUnh dtd " 


FI(7) 


1 1 



W(z-lr)t) 2 
Q (2n+5)/2 


dtdif 


FI(8) 


// 


w(z-int) 2 

Q (2n+7)/2 


dtdn 


w - n p4ttf2 (i-t I )“ 


Q = [ (z-irjt) 2 + r 2 ] 


(2.5.4) 


(2.5.5) 


(2.5.6) 


(2.5.7) 


(2.5.8) 


A computer program was written to integrate the FI(1) through FI (8) 
integrals numerically whenever they occur. The integration formula 
used was a 12th order Newton-Cotes closed type formula (21) . Even with 
this high order formula, a very fine mesh was required for small z and 
r < 1. In addition, the integration accuracy deteriorates with increas- 
ing n . However, despite these problems the integrals could be 
calculated to 3 digit accuracy for z as low as .2 for r < 1, and 
for n and p as high as 5. The only drawback in going to smaller 
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z and larger n , is that the required computer time tends to increase 
very rapidly. No integration problems were encountered when r > 1 , 
even for z^O. The stress components may be written as follows: 

N P 

a - l A(n) cos (n9) l C P H p r 11 SRR (2.5.9) 

rr n=0 p=0 n n 


N P 

O aa = I A(n) cos (n6) Y C p H p r n S06 (2.5.10) 

bo n n 

n=0 p=0 


N P 

T 1 A(n) sin (n0) £ C p H P r 11 SR0 (2.5.11) 

r0 n=0 p=0 n n 


N P 

T = 7 A(n) cos (n0) 7 C P H P r* 1 SRZ (2.5.12) 

rz n n 

n=0 p=0 


N P- 

= “ I A(n) sin (n0) £ C P H p r 11 S0Z (2.5.13) 

UZ n=0 p=0 n n 

where 

SRR » - Cl • FI(7) + C14 • FI (2) + C2 • FI(1) 

- C3 • FI(1) + C4 • FI(3) - C5 • FI(4) 

+ C6 • FI (5) + G7 • FI(6) (2.5.14) 

S00 = - Cl • FI (7) + C14 • FI(2) - C2 • FI(1) 

- C8 • FI(2) + C5 • FI(4) + C9 • FI(5) (2.5.15) 

SR6 = C2 • FI (1) - CIO • FI(2) - C5 • FI(4) 

+ Cll • FI (5) (2.5.16) 

SRZ = C12 • FI (7) - C13 * FI (8) - C15 * FI (2) 


+ C16 • FI(3) 


(2.5.17) 
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S0Z 


C12 • FI (7) - C15 * -FI (2) 



F(2n+1) 



/F 2 2a+1 [r(n+l)] 2 

Cl 

= 

2(8 2 -2)(2n+3)(2n+l) 

C2 

= 

2n(n-l)/r 2 

C3 

= 

2 (2n+l) 2 

C4 

= 

2 (2n+l) (2n+3) r 2 

C5 

= 

2(8 2 -l)zn(n-l) (2n+l)/r 2 

C6 

= 

2 (8 2 -l) z ( 2n+l) 2 (2n+3) 

C7 

= 

2(f3 2 -l)z(2n+l) (2n+3) (2n+5)r 

C8 

= 

2(2n+l) 

C9 


2 (8 2 -l) z ( 2n+l) (2n+3) 

CIO 

= 

2n(2n+l) 

Cll 

a 

2(|3 2 -l)zn(2n+l) (2n+3) 

C12 

= 

2(8 2 ~1) zn(2n+l) (2n+3)/r 

C13 

S3 

2(e 2 -l)z(2n+l) (2n+3) (2n+5)r 

C14 

=1 

2 (B 2 -2) (2n+l) 

C15 

= 

2(8 2 -l)zn(2n+l)/r 

C16 

ss 

2(8 z -l)z(2n+l) (2n+3)r 


(2.5.18) 
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The equations (2.5.9) through (2.5.13) together with the 
integrals given in equations (2.5.1) through (2.5.8) give the stresses 
near a circular crack in an infinite solid. Due to the problems with 
the numerical integration procedure for small z the special solution 
for z=0 and r < 1 will be presented in the next section. 

2.6 The Solution for z=0 and r < 1 

The stress components for the special case z=0 have the form 


rr 


z=0 


- 2(g z -2) 2 

3z 2 3r 2 


( 2 . 6 . 1 ) 


00 


z=0 


= - 2(0*-2) (2.6.2) 


■3z 


2 r 3r 2 


r 2 30 2 


r0 


2 3 2 (j> _ _2_ 3£ 

r 3r30 2 30 

z~0 r 


(2.6.3) 


The potential function itself becomes 


♦z-0 


- | y c p h p f n p+w do x 

n=0 JT n-“o n n J 


00 

/ 




/r 


(2.6.4) 


while -^| 
3z' 


z=0 


is given by 
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d 2 <t> 

dz 2 


z=0 


N 


I 

n=0 


cos (n8) 

/2 


P 

I 


p=0 




x 


/ W 5 " 5 V 5r) dE 

0 


(2.6.5) 


To determine the stresses on the plane z=0 , the integrals involved 
in these two expressions must be evaluated. 

The infinite integral in equation (2.6.4) is a Weber-Schafheitlin 
discontinuous integral. Abramowitz and Stegun (21) give the value for 
this integral as 


00 

/ 


J (at) J (bt) 

u V 


dt = 


b v r 

u+v-X+l 

3 

2 

2 X a v-X+l r(v+1) p . 

U-V+X+l 

2 


F f H+V-A+l y-v-X+1 . . b^ : 

o » O 9 vrX f 


for 0<b<a,X>-l and (y+v-X+1) > 0 , or 


a * 1 r 

fu+v-x+i 

l 2 


2 X b u-X+l r(u+1) r 

V-U+X+r 
2 . 


2 ’ 2 » » 2 


2-i 


for 0 < a < b , X>-1 and (p+v-X+1) > 0 

where the function (a,b;c;z) is the Gauss hyper geometric 

series (21). 



21 


The infinite integral of (2,6.4) becomes 


/ W 5t,) V 5r) f " 
0 ^ 


r 11 r (n-f%) 

ST n n+ ^ r(n+i) 


for 0 < r < ri , and 


n n+% y r(n+ s +%)r(s+%) [Vf 

Si r n+1 TT S=0 r < n+S+ 3/2) Si [yj 


for 0 < ri < r . 

(2.6.6) 


Substitution of equation (2.6.6) into equation (2.6.4) gives 


^ z =0 


y £9.g...M) y c p H p x 
n=0 Si p=0 n n 


/ T ' p+3 ' 5 { 

J 1/2 r”* 


y r (n+s+%) r (s+%) f n] 2s ] 

4 r/'r,4-o-4-,v. 1 i o j r f 0,1 


Si r n+1 tt s~ 0 r(n+s+3 ^ s 


f rp +3/2 1 

'l 1ST n 4 ^ r(n+i)J 


(2.6.7) 


Interchanging the integration and the summation in equation (2.6.7) and 
proceeding to evaluate the next integral gives 


4> 


z-0 


I c ° s (n6) I C P H P J rP +2 BT(n,p) 
n=0 ST p=0 j 

, W) (r n -rP"- 2 ) for p-n+1 # -1 

+ r(n+l ) \ 

j- r n log (r) ' for p-n+1 - -1^ 

(2.6.8) 
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where 


BT(n,p) = l 
s=0 


r (n-fs+%) r Cs-f%) 

r (n+s+3 2 )r (s+l) (n+p+2s+3) 


Equation (2.6.8) together with (2.6.9) give the value of <j> on the 
plane z»0. 

To find 3 2 <}>/3z 2 , the integrals in equation- (2.6.5) must 

Z— u 

be evaluated. A careful examination of the infinite integral in 
equation (2.6.5) shows that it does not exist. However, the form of 
the infinite integral can be changed by integrating the finite integral 
on ri by parts. The basic differentiation formula for Bessel functions 
can be used to show that 


In { 5 W Cn) } * W 5n) <2 - 6 ' 10> 

The integral on r| can be rewritten to give 

1 1 

/ i p+3<! W 5l,)d ' 1 ’ - / < T ' p ' ta+1 > * 

0 0 

[- rf n+ ^ J n ^(^n)dn 3 (2.6.11) 

Integrating (2.6.11) by parts gives 
1 

y n p+3 ^ - 

0* 


J k (0 } n p4% 

1 + ^ P +n+1 ^ J — J n -%^ r, > dT1 

0 


(2.6.12) 
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Substitution of equation (2.6.12) into the expression (2.6.5) gives 

„ y cos (n6) y C P H P x 

3z 2 n n=0 /l p=0 n n 

z=0 r 

1 00 

{ (p+n+1) / n p4% dn y* vC J n-% (?n) J n (5r)d5 

0 0 

00 

f vTJ n _ % (5) J n (Cr)d£ } (2.6.13) 

0 

As was expected, these integrals are also special cases of the Weber- 
Schafheitlin discontinuous integral and can be evaluated to give 

00 

/ vf J ,(£) J (r£)d£ =0 for 0 < r < 1 (2.6.14) 

n— Ar n 

0 

CO 

/ Sz J .(£ T 1 ) J (r£)d£ = 0 for 0 < r < n (2.6.15) 

j tl—*> n 

0 

00 

/ ^ W 5n) V rC>d? 

0 

Because the solution developed in 2.5 is perfectly acceptable for z=0 
and r > 1 the solution in this section will be developed only for 
the case r < 1 . With this restriction, 3 2 4>/3z 2 l Q may be written as 


. Ji n a ~ % Cr 2 -n 2 )~* 

r n /f 


for r > T} > 0 
(2.6.16) 
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9 2 4 > 


dz' 


2=0 


? C ° S ■■ ? c P H P X 

L n pr L n n n 

n=0 v2 p=0 


(p+n+1) / tv 

0 


4 . 

/ 


p+n 


r n /rr (r^-rf) 


2 _„ 2 \% 


dr) for r < 1 

(2.6.17) 


After performing the final integration, 9 2 <{)/9z 2 
9 2 <j> 


z=0 


takes the form 


9z“ 


2=0 


! 22S-M2. -I <p + " + 1 > C p h ' 

n=0 p=0 Al * i 


/rT 


for n+p = 0 


p n+p +2 

— r , o -i r P for n+p ^ 0 

l r 


for r < 1 


(2.6.18) 


The stress components for the plane z=0 can be calculated by 
substituting the equations (2.6.8) and (2.6.18) into the equations 
(2.6.1) through (2.6.3). 


rr 


N / QN P 

= y cos ■ (n - 6 . ) .. y c P H P r p X 
2=0 n=0 u p=0 n n 

r 

{ /n for n+p = 0 

H p for n+p 4- 0 

* n > 


+ < 


'"W) 1 "’' 2 - (p+MCp+D) 

for p-n+1 ^ -1 

- \(n+l)^ r 11 P 2 [n(n-l) log(r) + 2n - 1] 




(2.6.19) 


for p-n+1 = -1 


J J 
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? ^--. (n8) ! C P H P r P X 

L n IT L n n n 

n=0 p=0 


(p-n z +2) BT(n,p) - (0 2 -1 )tt 


vir for n+p = 0 


H y for n+p ^ 0 
n ' 


-LgfeW- [„(i- n)r ”-P - 2 + ^_ p . 2j 


r(n+X) (p-n+2) 


for p-n+1 ^ -1 


IT T (n+%) n-p-2 r ,, s , , „ . , , 

"TT n +l)~ r [n (1-n) log (r) +1] 


for p-n+1 = -1 

J 

( 2 . 6 . 20 ) 


T r0 „ * ? 2iS i Ml ! C n P H „ ” r ” j’ “<*•» 

2=0 n=0 p=0 . | 


- IT T (n+%) n r , i ^ n-p-2 , 

' f(n+l) (p-n+2) t( ”- 1)r -P- 11 f° r P--+ 1 * “I 


TSt+ffi n < n - 1 H 1 °8 M + U 


for p-n+1 = -1 

( 2 . 6 . 21 ) 


T = T Q =0 

rz - z8 _ 

2=0 z =0 


( 2 . 6 . 22 ) 


Equations (2,5.9) through (2.5.13) together with equations 


(2.6.19) through (2.6.22) constitute the complete circular crack 


solution 
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2.7 The Stress Intensity Factor 

The stress intensity factor is defined by 


K t = Limit /2tt6 o 

1 6 + 0 


(2.7.1) 


2=0 


where 6 is the distance from the crack tip and O is the normal 

r zz 

stress on the plane of the crack. Substitution of a into equation 

ZZ 


(2.7.1) yields 


N P 

K r = 2G(3 2 -1) l cos (n0) l E^ (2.7.2) 

n=0 p=0 

This method of calculating K^. was first used by Smith (13). 

2.8 Summary 

In Chapter II the problem of the circular crack in an infinite 
solid has been formulated and solved. Section 2,5 presents the stress 
components in the form of a double series. Some of the required 
integrations must be performed numerically and these are given in 
equations (2.5.1) through (2.5.8). Section 2.6 presents the stress 
components for the special case z=0 and r < 1. This special case 
had to be solved separately because of problems with the numerical 
integration procedure for small z . Section 2.7 presents a method for 
determining the stress intensity factor for an arbitrary normal pressure 
on the crack surface. 

The solution presented in this chapter can be used to calculate 
the stresses surrounding a circular crack in an infinite solid, when 
the crack surface is loaded with a normal pressure approximated by a 
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truncated Fourier-power series. The Fourier series contains terms up 
to cos (50) , while the power series includes terms up to r 4 . In 
addition, a direct method of calculating the stress intensity factor 
due to this loading has been presented. 



CHAPTER III 


THE HALF SPACE SOLUTION 
SOLUTION 2 


3.1 The General Problem 

In section 2.2 It was explained that solution 2 would he used 
to calculate the stress on the crack surface produced by freeing of 
the front and back surfaces during the iteration procedure. Figure 4 ■ 
is an illustration of how the front and back surfaces are arranged for 
this procedure. The following paragraphs describe how this portion of 
the procedure functions. 

After the stresses at the center of each rectangle have been 
computed using solution 1, solution 2 is used to calculate the stress 
on the crack surface produced by the freeing process. To accomplish 
this the crack surface is represented as a circular array of points 
and the stress at each of the array points is calculated by summing 
the stresses - produced by each individual rectangle of the front surface. 

The freeing stress which produces a stress a on the crack 
surface also produces stresses on the back surface. This is handled 
by considering each back surface rectangle one at a time and summing 
the stresses produced by freeing each of the front surface rectangles. 
Here again, the stress calculation is made for the center of each 
rectangle. The calculations proceed in an identical manner for iter- 
ation between the back surface and the crack. 
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To save computer time, the stress calculations between the 
front and back surfaces are made only after an iteration cycle between 
the crack and one surface has been completed. For example, if three 
iterations are to be performed between the front surface and the crack 
then the freeing stresses for each iteration are accumulated in a 
running sum and the back surface residual stress is calculated only 
once. 


3.2 The Solution for One Rectangle 

For any given rectangle the problem reduces to finding a 
solution for a semi-infinite solid when the rectangle on the boundary 
is subjected to uniform normal and shearing stresses. This basic 
problem has a rather long history and was once called the "problem of 
Boussinesq and Cerruti." A potential formulation for the problem is 
presented in Love's "A Treatise on the Mathematical Theory of Elas- 
ticity" (22). Later, in a paper (23) Love presents the solution for 
the case of constant normal stress. Much more recently Smith and Alavi 
presented the solution for the case of constant shearing stresses. 

Since the solution is rather involved and adds little to the under- 
standing of this work, the solution will not be reproduced here. For 
a more complete discussion of this solution and its application to 
similar problems in fracture mechanics see references (13) (14) and 
(15) . The formulas given by Love (23) , and by Smith and Alavi (15) 
were programed for the digital computer so that wherever the freeing 
process was required these subroutines could be used. 




CHAPTER IV 


THE COMPUTER PROGRAM 


4.1 The Program Logic 

This section gives a description of how the computer was 
programmed to use the two elastic solutions to solve the surface flaw 
problem. The essence of the material which will be given here was 
previously, presented in section 2.2. However, at that point it was 
intended as an introduction to the solution method rather than a 
description of how the solution was obtained. 

A block diagram of the computer program is shown in Figure 5. 
The number to the right of a block is the statement number in the main 
program where that block operation takes place. Although the block 
diagram does not give all of the details of the program operation, it 
does present the basic program logic. 

4.2 Program Checkout 

To insure that the program had no computational errors, a 
variety of test cases were run. For example, the circular crack 
solution was checked by comparing the results of this program with the 
closed form solution of Sneddon for the special case of a constant 
crack pressure. In addition, the program was checked against Alavi’s 
solution for a crack pressure described by three terms in the Fourier 
series. In both cases the comparison was accurate to at least three 
digits . 
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The overall results of the program were checked against the 
results of Smith (13) for the case of semi-circular crack, and against 
Alavi's results (14) for the cases D=.3 and D=.4 . The comparison 
with Smith is shown in Figure 6. The difference in results for 8/6^ 
near unity is about 5%. Most of this difference can be attributed to 
the fact that the front surface grids were not- the same. Smith used 
over 500 rectangles on the front surface while this study used 184 
rectangles for the first iteration and 62 rectangles for subsequent 
iterations. The close agreement between these two studies tends to 
indicate that the iteration procedure is not very sensitive to the 
surface grid used. The comparison with Alavi’s results was not plotted 
in Figure 6 because of the close agreement. For the case D=.3 the 
results were almost identical, and for D=.4 , Alavi’s results were 
about 2% higher than the results of this study. 

4.3 Surface Grids 

As mentioned above, two different front surface grids were 
used. In all cases the 184 rectangle grid was used for the first 
iteration, and the 62 rectangle grid was used for the following iter- 
ations. It was possible to use Sneddon's closed form solution to 
calculate the front surface stresses on the first iteration because 
the crack pressure is constant. Sneddon's closed form solution computes 
stresses at a much faster rate than the numerical method of Chapter II. 
For this reason, a finer grid could be used without excessive computer 
time. When Sneddon's closed form solution was used, the run time was 
about 200 seconds for this iteration. In addition, it is highly 
desirable to perform the first iteration as accurately as possible 
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since the stress intensity factor changes by the largest amount on the 
first iteration. For example. Smith found that for the semi-circular 
crack the stress intensity factor changed by 14% as a result of the 
first iteration and by less than 3% on the second iteration. The 184 
rectangle grid is shown in Figure 7, and the 62 rectangle grid is 
shown in Figure 8 and Figure 9, 

The back surface which has 32 rectangles is shown in Figure 10. 
The back surface is always located outside of the crack and never 
intersects the crack surface. For this reason the back surface is in 
a region of low stress gardients, and a fine mesh is not required. 

Smith and Alavi (15) calculated the stress intensity factor due to a 
free surface located outside of the crack surface using 64 rectangles 
on the free surface. To check the 32. rectangle grid some of the cases 
presented by Smith and Alavi were rerun. The difference between the 
results was virtually zero, 

4.4 The Fourier Series for Crack Pressure 

The crack pressure is approximated by a Fourier series in the 
0 direction and a power series in the r direction. See equations 
(2.3.3) and (2.4.1). As previously mentioned, the crack surface is 
represented as an array of points where the crack residual stress has 
been calculated at each of these points. The Fourier coefficients are 
determined by numerically integrating over the array’s 19 anular 
divisions for each of the 7 radial divisions, then a least-square curve 
fit is performed to obtain the power series coefficients of equation 
(2.4.1). Terms of up to cos (56) were included in the Fourier series 
and terms up to r 4 were used in the power series. 
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In order to realize the benefit of the 184 rectangle grid, and 
the closed form solution for the first iteration, . the stress intensity 
factor must be calculated with maximum accuracy. ' The equation for 
stress intensity factor (2.7.2) shows that this can be accomplished 
by increasing the number of terms in the crack pressure series. For 
this reason the Fourier series used terms up to cos (100) and the 
power series included terms up to r 5 for calculating the stress 
intensity factor due to the first iteration. This improved approxi- 
mation of the crack pressure was not used to compute the stresses for 
the next iteration, because the circular crack solution can only handle 
terms up to cos (50). This means that after the first iteration the 
crack pressure was approximated with a Fourier-power series containing 
terms up to cos (50) and r 4 . This less accurate approximation for 
the crack pressure is considered adequate because the second iteration . 
changes the results by less than 3%. 



CHAPTER V 


RESULTS 
PART I 


5.1 General Remarks 

For all of the results contained herein the iteration cycle 
was as follows: 

1. One iteration was performed on the front surface using the 
184 rectangle grid of Figure 7. The stress intensity factor 

_ was calculated using 11 terms in the Fourier series and 
6 terms in the power series. 

2. A second iteration was performed on the front surface using 
the 62 rectangle grid of Figure 8 and Figure 9. The 
pressure distribution was approximated with 6 terms in the 
Fourier series and 5 terms in the power series. The stress 
intensity factor was calculated using this same series. 

For the cases D = 0, .3, .5, .7 the iteration cycle was 
stopped at this point to give the solution for the circular 
crack partially embedded in a half space. 

3. For D = .2, .4, .6 the iteration cycle continued with 
two iterations between the back surface and the crack using 
the grid of Figure 10. The Fourier series contained 6 
terms and the power series had 5 terms. 


34 



35 


4. The iteration cycle continued with one additional iter- 
ation on each surface and then stopped. 

It is possible to approximate the semi-elliptical surface flaw 

I 

with a partially embedded circular crack. This can be accomplished 
several different ways, but in this work the comparison will be based 
on matching the curvature of the ellipse to that of the circular flaw 
for the same crack depth "A". This comparison is shown in Figure 11. 
This approximation will be used later to compare the results of this 
study with experimental data. 

5.2 The Results 

Figure 6 shows the stress intensity factor for the partially 
embedded circular crack. The calculated points through which the 
curves are drawn indicate an increasing tendency to wiggle with in- 
creasing D. This would be expected since is approximated by a 

truncated Fourier series and as D increases 0 decreases, and 

max 

therefore the lower order harmonics have less effect. The curves 

show the anticipated trend of decreasing stress intensity factor with 

decreasing crack depth. The trend of decreasing stress intensity 

factor with increasing 0/0 for higher D values is not unexpected. 

max 

Smith (24) noted the same trend for the semi-elliptical surface flaw. 

It would be dangerous to put too much faith in the results for ®/® max 
near unity because this region is one of high stress gradients due to 
the singularity at the crack tip. In addition, the surface grid in 
this region is relatively coarse. The comparison of the D=0 case 
for this study with the results of Smith (13) also tends to indicate 
an increased error in this region. However, the area of greatest 
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interest is near 6=0 , because this is the point at which the maximum 
stress intensity occurs for most crack geometries. 

Figure 12 through Figure 14 shows the effect of the back sur- 
face for D values of .2, .4, and .6 respectively. The results show 
the expected trend of increasing stress intensity factor as the back 
surface moves closer to the crack tip. 

Figure 15 is a plot of the back surface intensification factor 
plotted as a function of D and A/T . This factor is defined as 


M 


t 


*1 

(K^. at A/T = 0) 


for 0=0 (5.2.1) 


The real test of an analytical study is a comparison with 
experimental data. This comparison is presented in Figure 16 and 
Figure 17 . The data was taken from the work of Larson (25) . The 
material used in Larson's work was a brittle epoxy into which semi- 
elliptical surface flaws were placed. The comparison between the ex- 
perimental data and the theory is excellent. 



CONCLUSION TO PART I 


The problem of the circular crack partially embedded in a 
solid of finite thickness has been solved. The stress intensity factor 
has been determined for a variety of depths, D, and A/T ratios under 
uniform tension loading. During the course of the solution a numerical 
integration technique was devised to calculate the stress in the, 
neighborhood of a circular crack where the crack surface is loaded by 
a normal pressure P(r,8) which can be specified by a Fourier-power 
series in which terms up to cos (50) and r 4 are retained. In 
addition, the stresses on the plane of the crack have been determined 
in closed form. 

The stress intensity factor has been determined as a function 
of position around the crack front for the crack embedded a depth D, 
where D varied from 0 to .7 and for A/T ratios from 0 to .85. The 
results indicate that for deep cracks the back surface can raise the 
stress intensity factor by as much as 30%. These results are compared 
with experimental data for semi-elliptical surface flaw cracks. The 
comparison has been made by matching the curvature of an ellipse to 
that of the part-circular crack for a given crack depth A. The ex- 
perimental data scatter about the theoretical results and give a good 
comparison. 

Although all of the results presented here have been for the 
special case of Poisson's ratio equal to .25 the program is set up in 
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such a manner that other values of Poisson's ratio can be used. In 
addition, the program is not restricted to the special case of uniform 
tension. Any loading which can be approximated with the truncated 
Fourier-power series can be run with the program. 



PART II 


INTRODUCTION 

There are solutions to many two dimensional crack problems in 
linear elastic fracture mechanics (5) . These solutions generally in- 
volve the requirement that the loading causes the crack surfaces to 
move apart. In this portion of the thesis a solution for a two 
dimensional through crack is presented in which the requirement of 
crack surface separation has been relaxed. This type of fracture 
problem does not appear to have been discussed to any great extent in 
the literature. Burniston (26) has solved the partially closed 
Griffith crack for the case where the crack is closed at its center. 
However, a crack which is closed at one end is of more practical 
importance , since this could occur with a crack located in a bending 
field. A typical example of this situation is shown in Figure 18. 

Consider the specific problem of a crack located at the center 
of a pure bending' stress field. For this case the existing theory can 
be used to obtain an upper and lower bound for the stress intensity 
factor. The upper bound is found by loading the crack shown in 
Figure 18 with the linear load P' (y) * a i a Q ~ Ot x y '» where is the 

bending load intensity in psi/in. For this loading the stress intensity 
factor can be calculated using equation 33 on page 37 of reference (5) . 
This computation gives 

— ~2 ^i a o /lTa o 
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A lower bound for the stress intensity factor can be found by assuming 
that the crack of Figure 18 is open over one half of its length. The 
lower bound for the stress intensity factor is given by 


K, 


1 3 

2 ^ 2 


°4 a O 


The upper and lower bounds differ by a factor of 2 v'T . In Chapter 9 
this problem will be solved using the theory presented here and it will 
be determined that the actual stress intensity factor is about 80% of 
the upper bound. 

It is the objective of this portion of the thesis to determine 
the stress intensity factor for a partially closed Griffith crack 
subjected to a general load specified by a polynomial. This solution 
will be restricted to the case where the crack has only one open 


region. 



CHAPTER VI 


THE PROBLEM 

Consider an infinite two-dimensional, linear elastic, homo- 
geneous, isotropic material which has a through Griffith crack. The 
loading on the material is specified in the form of a polynomial such 
that the crack opens only once. The problem geometry is shown in 
Figure 18. 

The boundary conditions for this problem may be written as 

follows : 

1.0 =0 inside the open portion of the crack; 

XX A 

x=0 

U =0 outside the open portion of the crack; 

x=0 

T =0 all along x=0 since only symmetric 

^ x=0 loads will be considered. 

2. For large x and y all stresses must approach the 
loading stress P'(y). 

3. The stress at the closed end of the crack is finite. 

A classical approach to this type of problem is to remove the 
crack and solve the resulting elasticity problem. Then a second 
problem is formulated such that when these two problems are super- 
imposed the solution to the original problem is obtained. 
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Figure 19 illustrates how this technique is used to solve the 
problem of the partially closed crack. Problem A is the problem which 
results from removing the crack from the original problem. In the 
sketch of problem B it can be seen that the loading - P(y) is 
applied to the area where the crack is to be located. In problem A, 

+ P(y) is the stress at this same location so that superposition then 
gives the boundary conditions of the original problem. 

Problem A is a typical elasticity problem and can be solved by 
many different methods. Problem B is a mixed boundary value problem 
and its solution will be developed in the next chapter. 



CHAPTER VII 


THE MIXED BOUNDARY VALUE PROBLEM 


7.1 Boundary Conditions 

The boundary conditions for the mixed boundary value problem 
can be stated as follows: 


a. 


b. 


c. 


d. 


x=0 

- - P(y) 

for 

|y| 

< a 

(7.1.1a) 

x=0 

= 0 

for 

|y| 

> a 

(7.1.1b) 

x=0 

= 0 

for 

all 

x=0 

(7.1.1c) 

-> 0 

as x -*■ 00 




(7.1. Id) 


e. On an intuitive basis, the boundary condition at the closed 
end of the crack is 


= Limit /2Tr(y-a) 
y •+• a 


a 

XX 


= 0 

x=0 


(7.1.1e) 


This supposition is easily verified by considering the alternate 
possibilities: Suppose that Rj. < 0, then calculation of the crack 

opening displacement shows that the crack surface has deflected through 
itself, which is impossible. Suppose that > 0, then from the 
definition of stress Intensity factor the stress cr > 0 at x=0 and 
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y=a, which is also impossible because the physical crack extends past 
y=a and cannot support a tension stress. The only remaining alternative 
is that Kj. = 0. The problem geometry is illustrated in Figure 20. 

7.2 The Potential Formulation 


For the case of plain strain the Navier equations take the 

form 



\7 2 TJ 4- ^ 

v u + (1-2V) 

3_ 

3x 

{ 1 + 

3V 'l 
3y J 

► o 0' 



V 2 V + ^ 

V V + (1-2V) 

3_ 

3y 


3V -1 

3y J 

\- 0 

(7.2.1) 

where 

U = displacement 

in the 

X 

direction 



V = displacement in the y direction 
v = Poisson’s ratio. 


For the class of problems in which the shearing stresses vanish at all 
points on the plane x=0 the following potential formulation satisfies 
Navier* s equations (7.2.1). 

u - - 8 2 || + (g 2 -i) x 

3x 3x 2 

v = |i+ ( 82 - 1 ) x Sfi (7.2.2) 

3y 3x 2 

o2 _ 2(l-v) 

P (l-2\» 

Provided that cf) satisfies V 2 <j> = 0 

The above displacements U and V have been nondimensionalized with 
respect to the crack radius "a". 
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Using Hooke’s law the stress components in terms of the poten- 
tial function <j> are found to be: 


CT - - 2(3 2 -l) ^ + 2(e 2 -l> x ^ 
XX 3x 2 3x 3 


a = - 2 ($ 2 -l) ^ - 2(g 2 -l) x ^ 

™ 3x 2 3x 3 


(7.2.3) 



2(f5 2 ~l) X 

3x 2 3y 


These stresses have been nondimensionalized through division by G , 
the shear modulus. 

The boundary conditions (7.1.1a&b) now take the form 


= P(y) 

3x 2 2(3 2 -l) 


|y| s a 


(7.2.4a) 


3 £ 

3x 


0 


y| > a (7.2.4b) 


The boundary condition (7.1.1c) is automatically satisfied by the 
potential formulation. 


7.3 The Potential Function 

To satisfy the mixed boundary conditions of equations (2.2.4a 
and b) a potential function will be assumed in the form 
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4> 


1 

2 (3 2- l) 


J f ^ - cos (?y) e ^ x 
0 


dc + ; 


2(e 2 -D 


X 



sin (5y) e~ Cx d£ 


(7.3.1) 


Potential functions in this form have been used by Sneddon (27) to 
solve related problems in fracture mechanics. From the form of cj> it 
is clear that Laplace's equation is satisfied, and that for large x 
values a., tends toward zero. 

Substitution of the above <j> function into the boundary con- 
ditions (2.2.4a and b) gives the following set of integral equations: 


CO 

/ 


P(y) = / Z f(5) cos (£y) d£ 

"0 


CO 

+ J £ g(D Sin (£y) d£ , |y| £ a 


(7.3.2) 


CO 


f(5) cos (£y) d£ 


CO 

/ 


+ / g(£) sin (£y) d£ , jy| > a . (7.3.3) 


Now if P(y) is divided into a symmetric and an, anti-symmetric function, 
the above integral equations may be rewritten as a pair of dual integral 
equations . 
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The symmetric portion gives the following set of dual integral 
equations : 


CO 

P s (y) = J K f(£) cos (£y) d£ , I y I < a (7.3.4) 


oo 

0 = f f(£) COS c^y) d£ , |y| > a (7.3.5) 


where P (y) is the symmetric portion of P(y). 
s 

The anti-symmetric portion gives the pair of dual integral 
equations : 


OO 

? a (y) “ J K g(?) sin (£y) d£ , Jyj < a (7.3.6) 


0 = J g(?) sin (£y) d£ , |y| > a 

0 


(7.3.7) 


where P (y) is the anti-symmetric portion of P(y). 

The solution to the above dual integral equations, due to 
Busbridge, is given by Sneddon in (18). 

The relationships which give f(£) and g(£) are 


2 r r V p) 

f (5) - f J r,J 0 (&i)dn J 7 ===== dp 

n n y r) ~ P 


(7.3.8) 


P * (P) 


g(?)=f / JjttXJdX f 


dp 


(7.3.9) 



With the functions f(£) and g(£) known, the potential 


function (J) is completely determined except that' "a", 
portion crack length is still unknown. To determine " 
condition 1^ = 0 at the closed end of the crack (2.1 
used. However, this will not be done until after the 
component has been computed. 


the open- 
s'', the boundary 
•le) will be 

cr stress 
xx 



CHAPTER VIII 


THE STRESS AND DISPLACEMENT 


8.1 The Stress and Displacement Functions 

To determine the stress and displacement on the plane x=0 
the functions f(£) and g(£)> (7*3.8) and (7.3.9), are substituted 
into the potential function $ and then the formulas (7.2.2) and 
(7.2.3) are used to calculate the stress and displacement. This 
computation gives 


XX 


x=0 


oo a 


-/{*/ 


n j q (£ti) dp x 


0 v 0 


/ 


r P <P> 


A] 2 - 


} 


dp > £ cos (£y) d£ 


oo a 


/{' / J ' (w " 1 dp l x 


0 v 0 


£ sin (£y) 


for [ y | > a (8.1.1) 


and 


U 




x=0 2(p 2 -l) 


/{» h 

0 1 0 


J n (?n) dn x 


'i 

/ 


v<» 

/n 2 - i 


} 


dp >• cos (£y) d£ 
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3 2 

2(3 2 -!) 



Jj(?X) dX 


x 


/ 


P P a <P) I 

- • dp L sin (£y) d5 

/x 2 - p 2 J 

for |y| 


< a 


( 8 . 1 . 2 ) 


Before the above integrations can be completed the form of the loading 
function P(y) must be known. To make the integrations as simple as 
possible the loading function will be assumed in the form of a poly- 
nomial . 


8.2 The Loading Function P(y) 

The loading must be specified as some function in the x,y 
coordinate . system of Figure 18, However, the loading is known only in 
terms of some material reference frame. For this problem the material 
reference frame is located at the center fo the physical crack and is 
designated x* , y* in Figure 18. In terms of the x* , y* coordinate 
system the loading function may be written as 

P(y') = a Q + otjy' + a 2 y’ + . . . + a^y’ 

K 

p(y») = £ cty’^ (8.2.1) 

n=0 n 

From Figure 18 it is clear that 

y’ » y - (a 0 “ a > 

With this relationship the loading function may be written in the 


form 
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K 


P(y) = l 


m=0 


3 y 

m 


m 


( 8 . 2 . 2 ) 


where 


3 m = I a n ('§) (a - a G ) 


n-m 


n=m 


n f. 


(8.2.3) 


and (m) are the binomial coefficients. 

Now the loading function may be broken into an even and an odd 
function to give 


P s (y) + p a (y) - f s 2m y 2n + 1 ^y 2 " 1 

m=0 m=l 


(8.2.3) 


where H+l is the number of even terms in the loading function and J 
is the number of odd terms. 

8.3 Integrated Form of the Stress and 
Displacements 

Now that the polynomial form of P(y) has been established the 
expression for stress and displacement can be determined. Denoting the 
second term of -equation (8.1.1) by "I" and then rewriting it gives 


y r r p p * (p) dp 
i = -?/ dx / 


0 0 0 


oo 

J 5 Sin (?y) J x (a> d£ (8.3.1) 


Making use of the relationship 




- j x (a> 5 


(8.3.2) 


The infinite integral may be rewritten as 
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d . /V 

-}/<*/ 


p P (p) dp 
a 




oo 

{•i/ 

^ A 


J Q (£X) sin (£y) d? 

(8.3.3) 


} 


Watson- (20) gives the value for this infinite integral, which is a 
special case 'of the Weber-Schafheitlin integral, as 


/ 


J 0 (5A> sin (£y) d£ 



for X < y (8.3.4) 


Substitution of this result into equation (8.3,3) gives 


I 


2 f X d X 

* { < y 2 - x 2 )^ 2 



p P (p) dp 

A 2 - p 2 


From equation (8.2.4) P (p) 

d 


can be written 


(8.3.5) 


■ P a ( P> 


I 


m=l 


g 2m-l 


2m-l 

P 


(8.3.6) 


Substitution of P (p) 
a 


into (8.3.5) gives 


I 



X d X 

(X 2 - P 2 ) 32 


J- 


I 3 

m^l 


2m- 1 


/ 



dp 

(8.3.7) 


The last integral in (8.3.7) has the value 


/ 



/rF r(m4%) ,2m 
2 T(m+1) A 


2m > 1 (8.3.8) 
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Substitution of this result into the expression (8.3.7) results in 
the term 


I 


1 y p T (m4%) 

w »-i 2 - 1 r(n+l) 


/ 


, 2m+l 

— — dA 

(y 2 -* 2 )’* 


(8.3.9) 


Using a similar procedure on the remaining terms in equations 
(8.1.1) and (8.1.2), it follows that the stress on the plane x=0 is 
given by 


a 

xx 


x=0 


a 

XX 


x-0 


B f r -^ 1 

° l 


H 


+ I 3 


r(m-H?) 


/it m=l 2rar ( m+1 > J (y 2 -a 2 ) 3/2 


f* 


2m+l 


dn 


. 2m+l 


. _ 1 _ r r(nri-%) r A dA 

/rr m=l 2m_1 r(m+1) ' (y 2 -A 2 ) 3 ' 


for yj > a 


P(y) 


for jy| < a 


(8.3.10) 


and the displacement on the plane x=0 


is given by the expression 


U 


x=0 


= a-v) 
/if 


f ? r(nrHf) f x 2m+1 

l Jo S2 “ r <“ +1 > 7 ^ 


+ m l / J=f d * } £or hi 5 3 


-u 


0 for y > a 


(8.3.10) 


lx=0 
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The following three integration formulas will be used to 
obtain the final expression for stress and displacement: 


/ 


2m+l 



/ 2 . 2 
/a -y 



k=0 


(a 2 -y 2 ) m - k 

[2(m-k)+l] 



for m > 0 


(8.3.12) 


y 


/ 


2m-l 
x 

r~i 2 

vx -y 


dx 



m-1 


l 

k=0 


, / 2 2Nin-k-l 

f™" 1 ) (a -y ) 

' k t [2(m-k)-l] 


2k+l 

y 

(8.8.13) 


for m > 1 


/ 


2m+l 
x 

(y 2 -x 2 ) % 


a 2(m+l) (2n+1) m (g) H) k 

y 2 y io t 2 <”- k > +1 > 

y 215 ' 1 - y 2 ”] (8.3.14) 

for m > 0 


Now the stress on the plane x=0 can be written as 


xx 


x=0 • y 2 /if ^y 2 -a 2 


K 


.3 /iF f g r (m4%) 2m- 1 

r ~ + y i e 2m r(m+l) a 

a* m=l 


+ y g r. M g). a 2m_1 1 _ B + _L f g x 

il ^ r <“ +1 > J B ° /F m-1 


r(m+%) 

r(m+l) 


(2m+l) 


? ( 11 ^ 

k £ 0 [2(m-k)+l] 


f / 2 2\ti-*k+% 2k 2m 

[ (y -a ) y -y 


i V ft r(nri») (2m+l) ? (k) < ~ 1) 

/jf m =! 2m- 1 T(m+1) y k f 0 [2(m-k)+l] 
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2 v m-k4% 2k 2m 


[(y 2 -a z r - y - y 


for y > a 


and 


a 


xx 


= P(y) for |y < a 


x=0 


The displacements take the form 


U 


x=0 


(1-V) 

/if 



r(nH%) £ 

r(m+1) k =Q 


(£) (a 2 -y 2 )“~y k 

I2(m-k)+l] 


and 


+ 


J 


l B 

m=l 


r(nH%) 

2m— 1 r (m+1) 



0 


m -1 ( m - 1 ) ( a ^ y 2 ) 1 *-^ 1 

k=0 [2(m-k)-l] 

for | y j < a 

for ) y j > a 


2k+l 

y 


} 


8.4 The Open Crack Length "a" 

The only remaining unknown is the open crack length "a". To 
determine "a" the last boundary condition, = 0 at the closed end 
of the crack (7.1,le) will be used. This condition may be written as 


Limit /2ir(y-a) 
y -*■ a 


O 

xx 


x=0 


0 


(8.4.1) 


To simplify the equations the stress will be written in the following 
symbolic form: 




XX 


x=0 


a ” 

/* y 2 /y+a. 
/y-a 


H(y) + G(y) 


(8.4.2) 
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where 


H(y) = 3 0 y 3 + y I 3 


r (m4%) 2m-l 
2m F (m+1) 3 


+ ^ ^2m-l 

m=l 


T (m+%) - 2m- 1 

r(m+l) 3 


(8.4.3) 


G(y) = 


ft + — V ft 


rcrt) (2m+1 . ? «) ( ~ 1)k 

r(m+l> Um+i) k ; 0 [ 2 (m-k)+l] 


2 2 N m-k-*% 2k-l 2m"l . 1 . v 0 v 

L (y " a } y - y J + - I fWi x 

/ 7T m=l 


T (m+%) (2m+l) y (k) (-1) 

T (m+1) y ,“ n [2(m-k)+l] 


[ ( y a -a 2 )" lkH% y 2k_1 - y 2m ] 


(8.4.4) 


With the substitution of the symbolic form for O into the 

xx 

boundary condition (8.4.1) the following requirement results 


Limit /2Tr(y-a) 
y -*■ a 


H(y) + G(y)f = 0 (8.4.5) 


Careful examination of equations (8.4.3) and (8.4.4) shows that both 
H(y) and G(y) are finite at y=a. Therefore, to satisfy equation 
(8.4.5) H(y) must be factorable into the form 


H(y) = (y-a) Q(y) 


If H(y) is factorable into the form shown in equation (8.4.6) then 
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H(a) 


Y o r(m-^) 

Jo 2m r(nl+1) 


2m 


J 

+ I 

m=l 


r(nri%) 

2m-l T (m+1) 


2m-l 

a 


0 


(8.4.7) 


and the open crack length "a" is found to be one of the roots of the 
equation H(a) = 0. 


8.5 Admissible Roots 

Some of the roots of H(a) = 0 are not admissible as values 
for "a". For example, complex roots have no physical meaning. Refer- 
ring to Figure 18, the following statements can be made: 

1. One real root on [0 , a Q ] and, or one real root on 

[0 , -a Q ] is allowable. This restriction must be imposed 
because the initial formulation was made assuming that 
the crack opened only once. 

■ 2. Roots where |a| > | a [ have no physical interpretation. 
3. Complex roots have no physical- interpretation. 


The solution to the partially closed crack is now complete, 
however, one additional point will be investigated before proceeding 
to an example. 


8.6 Slope at the Closed End of the, Crack 

For the Penny shaped crack, Emery and Smith (27) '“found that 
when Kj. = 0 the slope of the crack opening shape was zero at the 
crack tip. Similar results can be expected here. This can be written 
in equation form as 
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9U 

3y 


x=0 

y~ a 


0 


( 8 . 6 . 1 ) 


In this section it will be shown that this conjecture is true. 

The first step is to take the derivative of U with respect 
to y . This computation gives 


3D 

9y 


x=0 



R(y) + /a 2 -y 2 


( 8 . 6 . 2 ) 


where 


R(y) 


TT / TH \ / 2. 2 . \Hl“k 2k 

(i-v) ? o rwi ? It) Ca -y } y 

/T m=0 2 ” r( " +l) k=0 U(m-k)+l] 


(1-V) 

/F 


J 

I 

m=l 


6 


F (m+%) 
2m-l T(m+1) 


m-1 

I 

k=0 


rs 1 ) 


( 2 2\Bl-k-l 2k+l 
(a -y ) y 


[2(m-k)+l] 

(8.6.3) 


The next step is to examine the behavior of (8.6.2) as y approaches 
"a". Since the second term of (8.6.2) goes to zero this leaves 


3U 

9y 


x=0 

y=a 


Limit 

y a 


• -..I Limit { R(y) i 

.y^?\ y - a 1 . J 


(8.6.4) 


It can be demonstrated that 


Limit R(y) = 
y -> a /F 


f H 

I 

l m=C 


„ F (m4%) 2m 

e 2 m rSS) a 


+ I S 2m-1 

m=l 


F (m+%) . 2m- 

T(imfl) a 


*} 


(8.6.5) 
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However, the term in brackets Is H(a) which is known to be zero. By 
slimming out R(y) in equation (8.4.3) it can be shown that R(y) is a 
polynomial in y . Now if the polynomial R(y) is equal to zero at 
y = a , then R(y) is factorable into the form R(y) = (y-a)B(y) . 
Hence, the Limit (8.6.4) can be rewritten and evaluated to give 


3U 

3y 


x=0 

y=a 


Limit (1 V) y &~3). B(y) 
y -*■ a SiT /a 2 -y 2 


0 


( 8 . 6 . 6 ) 


This shows that for a loading function P(y) in the form of a polynomial, 
if Kj: = 0 at y = a , then 3U/3y =0 at y = a , 



CHAPTER IX 


AN EXAMPLE PROBLEM 

9.1 The Problem of Pure Bending 

Consider a bending problem where the applied load, P'(y), of 
Figure 18, has the form 

P'(y) ** a Q + a x y (9.1.1) 

For this applied load, the loading function P(y) of section 8.2 is 
given by 

p(y) = e o + 3 x y (9.1.2) 

where, from (3.2.3) 

£ = a + a. (a - a ) (9.1.3) 

o O 1 o 

= a x (9.1.4) 


The open crack length "a" is found by substitution into (9.4.7) 
which gives 


a 


2 

3 


(a, a - a ) 

A o o 


a 


i 


(9.1.5) 


For the special case of pure bending. 


a - 0 , and then a = 2/3 a . 
o o 
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9.2 The Stress and Displacements for Pure Bending 

To find the stress and displacements on the plane x = 0 
requires substitution of 3 q and into equations (8.3.15) and 

(3.3.16). This substitution gives 


a 

xx 


x=0 


M 


0 yCW 


- 1) 




-z 2|y| ) 

(9,2,1) 


U 


x=0 


<* 0 + 2 } 


(9.2.2) 


Great care must be taken to use the correct sign for the terms which 
result from the anti-symmetric loading. To calculate the influence of 
an anti symmetric term for y < 0 , y is replaced with |yj and then 
each anti-symmetric term is multiplied by minus one. Equations (9.2.1) 
and (9.2.2) have been setup in this manner. 

Figure 21 is a plot of the stress on the plane x^O. This is 
the stress for the mixed boundary value problem of Chapter VII, which 

is illustrated in Figure 20. 

* • 

Figure 23 shows the stress which results from the complete 
problem as shown in Figure 18. This is the final result for the 
original problem posed in Chapter VI. 


9.3 The Stress Intensity Factor 

Substitution of the stress (9.2.1) into equation (2.7.1) for 
stress intensity factor gives 


K = Limit ) 
y -> a 


r/2ie o 


l ^ 


± B, 


m y 2 
/y+a 


1 


(9,3.1) - 
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for a = 0 , then 
o ’ 



(9.3.2) 


Then at the closed end of the crack K^. = 0 and at the open end of 
the crack 





(9.3.3) 


The stress intensity factor given in equation (9.3.3) is about 80% of 
the upper bound, which was calculated in the introduction to Part II. 



CONCLUSION TO PART II 


A formulation for the partially closed Griffith crack has been 
presented. The resulting mixed boundary value problem has been solved. 
The following results have been obtained: 

1. A closed form solution has been obtained for the stress 
and displacement on the plane x=0. In addition, a 
criterion for determining the open crack length "a" has 
been presented. 

2. These results show that for the closed end of the crack 
the requirement that the stress intensity factor be zero 
aiso requires that the slope of the crack opening shape 
be zero at that point . 

3. To demonstrate the use of these results the problem of 
pure bending which closed one end of the crack has been 
solved. 

The general form of the results which have been presented here make it 
possible to solve an entire class of fracture problems where the crack 
is forced closed at one or both ends. 
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CONCLUSION TO THE THESIS 


The surface crack has often been found to cause the failure of 
structural components; however, as previously mentioned, only limited 
analytical work has been done on this problem. This is the first 
analytical study of a surface crack in a solid of finite thickness. 

The results of this study can be directly applied to the fracture 
analysis for brittle materials. This is accomplished by calculating 
the stress intensity factor from the results presented in Figure 6 
and Figure 15. - 

The design engineer attempting to design against a fracture 
failure is likely to select a material which is very .ductile. In that 
case plasticity effects would be significant and the results of this 
study would not be applicable. However, the results of this study 
together with an experimental program would make it possible to 
determine the plasticity effects. The results of the experimental 
program could take the form of a plasticity correction factor which 
would be applied to the elastic stress intensity factor of this study. 

For two dimensional cracks the size of the plastic zone around 
the crack tip has been estimated using the elastic analysis. This 
has been accomplished by determining the region in which the elastic 
analysis indicates that yielding has occurred. In a similar manner, 
the elastic analysis presented here could be used to approximate the 
yield zone for a surface crack in a finite solid. 
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Based on the success of this study it would seem worth while 
to use this same iteration procedure with the elliptical crack solution 
(11). In this way, the semi-elliptical crack could be studied directly. 

The solution which has been presented for the partially closed 
Griffith crack extends the theory of linear elastic fracture mechanics 
to include a new class of problems. The solution here makes it pos- 
sible to solve the class of problems where the loading closes a portion 
of the crack. The problem of pure bending which closes one end of the 
crack is the most common example of this situation. Although the 
problem of pure bending is the only example which has been worked out 
here, many different types of problems can be solved with the solution 
presented. The solution is general because the loading on the crack 
surface is represented in the form of a polynomial. For example, a 
thermal stress problem could be solved by approximating the crack 
surface stresses due to the thermal loading with a polynomial. Then 
this polynomial would be used as the loading function in the analysis. 

The general solutions presented in this thesis represent a 
significant advance in the theory of fracture mechanics. The results 
presented here are directly applicable to a wide class of problems, 
which up to this point have had no analytical solutions.- 
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Figure 2. — The Circular Surface Flaw 
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Figure 3. — The Circular Crack in an Infinite Solid 
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Figure 4. — Setup For the Iteration 
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Figure 5. — Block Diagram 
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Figure 6. — The Stress Intensity Factor for the Partially Embedded Circular Crack 
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BLOCK RECTANGLE SIZE 



Length yp : 

x Length zp 

A 

.05 

X 

.05 

B 

.1 

X 

.05 

C 

.1 

X 

.075 

D 

.1 

X 

.1 

E 

.2 

X 

.2 

F 

.5- 

X 

.5 

G 

1. 

X 

1 . 


Figure 7. — Front Surface Grid (184 Rectangles) 
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— Front Surface Grid ( 




— Back Surface Grid 


76 



A/2C « 1/2 /MD 


D 

A/2C 

0 . 

.5 

.3 

.42 

.4 

.39 

.5 

.35 

. 6 

.32 

.7 

.27 

.8 

.22 


Figure 11. --Matching Curvature of the Semi-Ellipse 
to the Part Circular Crack 
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Figure 12. — Back Surface Effect for D = .2 
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EXPERIMENTAL DATA FOR A/2C = .35 
COMPARED TO THEORY FOR D = .6 
(A/2C = .32 MATCHING CURVATURE) 
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Figure 20, --The Mixed Boundary Value Problem 
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NOMENCLATURE FOR PART I 

Stress Components 

Displacement Components 

Crack Pressure 
Crack Depth 
Thickness 
Shear Modulus 
Poisson Ratio 
Lame'' Constants 
2(1 - V)/(1 - 2v) 

Potential Function 
Fourier Coefficients 

Power Series Coefficients 

Bessel Function of the first kind of order v 

r( )/r< > 

Stress Intensity Factor 

Depth of Circular Crack 
Circular Cylindrical Coordinates 
Rectangular Cartesian Coordinates 
Distance from the Crack Tip 
Crack Radius 
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NOMENCLATURE FOR PART II 

Stress Components 

x Displacement 
y Displacement 

Length of Open Portion of the Crack 
Physical Crack Length 

Load on the Infinite Space 
Load on the Crack Surface 
Poisson Ratio 
2(1 - V)/(l - 2v) 

Shear Modulus 
Potential Function 
Symmetric Portion of P(y) 

Anti-symmetric Portion of P(y) 

Loading coefficients for the Polynomial P' (y) 

Loading Coefficients for the Polynomial P(y) 

Bessel Function of the first kind of order V 

Stress Intensity Factor 

Two Dimensional Rectangular Cartesian Coordinates 
Distance from the Crack Tip 
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